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Abstract 

Given the n x n matrix polynomial P(x) = YI^—q PiX 1 , we consider the associ- 
ated polynomial eigenvalue problem. This problem, viewed in terms of comput- 
ing the roots of the scalar polynomial detP(x), is treated in polynomial form 
rather than in matrix form by means of the Ehrlich-Aberth iteration. The main 
computational issues are discussed, namely, the choice of the starting approx- 
imations needed to start the Ehrlich-Aberth iteration, the computation of the 
Newton correction, the halting criterion, and the treatment of eigenvalues at 
infinity. We arrive at an effective implementation which provides more accurate 
approximations to the eigenvalues with respect to the methods based on the QZ 
algorithm. The case of polynomials having special structures, like palindromic, 
Hamiltonian, symplectic, etc., where the eigenvalues have special symmetries in 
the complex plane, is considered. A general way to adapt the Ehrlich-Aberth 
iteration to structured matrix polynomial is introduced. Numerical experiments 
which confirm the effectiveness of this approach are reported. 
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1. Introduction 

Given two positive integers k, n and matrices Pj € C" xn , j = 0, . . . , k 
consider the matrix polynomial 

p{ X ) = Y,p^ a) 

3=0 
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where Pk ^ 0, so that P(x) has degree k, and define the scalar polynomial 
p(x) = det P(x) of degree N < nk. Assume that P(x) is regular, that is p(x) is 
not identically zero. 

Under such an assumption, the polynomial eigenvalue problem associated 
with P(x) consists in computing the roots of the polynomial p(x) which are 
called the eigenvalues of the matrix polynomial P(x). Observe that, if Pk has 
not full rank, then N < nk. In this case, it is convenient to introduce nk — N 
eigenvalues at infinity and say that the matrix polynomial P(x) has nk eigen- 
values including the nk — N eigenvalues at infinity. 

Our interest is addressed to the design and analysis of efficient algorithms 
for the polynomial eigenvalue problem based on the Ehrlich-Aberth iteration 
0,0,0. 

Recently, much literature has been addressed to the polynomial eigenvalue 
problem (PEP). For the numerical solution of PEPs, fast and numerically stable 
methods are sought. Several algorithms have been introduced based on the 
technique of linearization where the polynomial problem is replaced by a linear 
pencil with larger size and the customary methods for the generalised eigenvalue 
problem are applied. For more details, see for instance [la 1221 l2a |47 1 an d the 
references therein. 

Specific attention concerns structured problems where the matrix coefficients 
Pj have some additional property which is reflected on structural properties of 
the roots. For instance, in the case of T-palindromic polynomials [23|,|43|, where 
Pj = PfrLj € C" x ™, the roots are encountered in pairs (x, l/x). In general, we 
may consider the case of structures where the roots can be grouped in pair s as 
(x, f(x)), where f(x) is any analytic function such that f(x) = f^ 1 (x) [l4j. In 
this case the goal is to design algorithms which take advantage of this additional 
information about the eigenvalues and deliver approximations to the eigenvalues 
which respect these symmetries independently of the rounding errors. 

The Ehrlich-Aberth iteration (EAI) was historically first mentioned in 0] 
and afterwards independently rediscovered many times. It is one of the many 
simultaneous iteration techniques available in the literature for the numerical 
approximation of polynomial roots 29|, 4(J ■ In 0, 0] the EAI has been combined 
with various techniques like the Rouche theorem, the Newton polygon technique, 
and the Gerschgorin inclusion theorems for arriving at efficient and robust soft- 
ware implementations. The package Polzeros, designed in [3|, provides a robust 
and reliable tool for approximating roots of polynomial in floating point arith- 
metic. The package MPSolve designed in Q provides certified approximations 
to any desired number of digits of the roots of any polynomial. 

The EAI has been used in [f| to solve the generalised tridiagonal eigenvalue 
problem where the software provides effective accelerations in terms of CPU 
time. It has been used in [4l| for quadratic hyperbolic tridiagonal eigenvalue 
problems. 

In this paper we present an adaptation of the Ehrlich-Aberth method for the 
numerical solution of PEPs. The main computational issues that we analyze 
are the choice of the starting approximations, the computation of the Newton 
correction, the halting criterion, the design of a posteriori error bounds, and the 
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management of the problematic (multiple) eigenvalues at zero and at infinity. 

Concerning the choice of the starting approximations, we propose a gener- 
alization to matrix polynomials of a technique introduced in [3j for scalar poly- 
nomials based on the Rouche theorem. In fact, we rely on the generalization to 
matrix polynomials of the Rouche theorem given in |34j | and provide a way to 
determine an annulus in the complex plane which contains all the eigenvalues. 

The Newton correction is computed by means of Jacobi's formula for the 
differential of the determinant of any square matrix in terms of the trace of the 
matrix P(z)^ 1 P' (z). The halting condition is given in terms of the condition 
number of P(z). Eigenvalues at zero and at infinity can be removed, to a certain 
extent, by using the specific features of the EAI relying on the information 
provided by the singular values of Pq and Pk- A posteriori error bounds are 
given by constructing a set of inclusion disks relying on the Gerschgorin theorem 
and adapting the results of 45] to the case of matrix polynomials. 

The computational analysis of this method shows that the number of arith- 
metic operations (ops) is 0(k 2 n 3 + kn 4 ). In the case where the degree k is large 
with respect to the square root of the matrix size n, this complexity bound 
compares favourably with the bound 0(k 3 n 3 ) of the customary matrix-based 
algorithms like the QZ applied to a linearization. Cases of this kind can be 
encountered, for instance, in the truncation of matrix power series (49j . 

The EAI does not compute the eigenvectors, which are sometimes needed as 
well. Nevertheless, after a good approximation of the eigenvalues is obtained, 
other methods, e.g. the SVD or the inverse iteration, can be used to compute 
the eigenvectors without increasing the complexity of the algorithm. We were 
able to compute eigenvectors with high accuracy using the eigenvalues given by 
the EAI. 

We consider the case of polynomials endowed with specific properties like 
palindromic, T-palindromic, Hamiltonian, symplectic polynomials, whose eigen- 
values have special symmetries in the complex plane. We propose a unifying 
treatment of this class of structured polynomials and show how the EAI can be 
adapted to deal with these classes in a very effective way. In fact, our variant 
of the EAI enables one to compute only a subset of eigenvalues and to recover 
the remaining part of the spectrum by means of the symmetries satisfied by the 
eigenvalues. By exploiting the structure of the problem, this approach leads 
to a saving on the number of operations and provides algorithms which yield 
numerical approximations fulfilling the symmetry properties. 

We conclude our discussion by presenting the results of several numerical 
experiments performed in order to test the effectiveness of our approach in terms 
of speed and of accuracy. We have compared the Ehrlich-Aberth iteration with 
the Matlab functions polyeig and quadeig[16j. In the structured case, we 
have also considered, when available, other structured methods, say, the URV 
algorithm by Schroder . 

We show that the EAI is much faster than the available techniques in the case 
where the degree is larger with respect to the size of the matrices. Moreover, for 
the test problems NLEVP of [2| , it turns out that the accuracy of the computed 
approximations is generally better than the accuracy obtained with the available 
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algorithms. 

The paper is organised as follows. In Section [2] we recall the Ehrlich-Aberth 
method and discuss the main computational issues encountered in its imple- 
mentation. In Section [3] we consider the case of "structured polynomials" , i.e. 
the case where the matrix coefficients have some special properties. Section 2] 
reports the results of the numerical experiments. 



2. The Ehrlich-Aberth method for matrix polynomials 

Given a vector £ <C N of initial guesses for the TV roots of the polynomial 
p(x), the EAI provides the sequence of simultaneous approximations given 
by 



jV(x) 



3 l-M{yf)My^ P'W 

( 2 ) 

where Af(x) is the Newton correction. It is easy to check that the jth update 
in ([2]) is nothing else but the Newton iteration applied to the rational function 

P( x )/YieLi e^j( x ~ Vi )■> so * n at the EAI provides a way to implement the 
implicit deflation of the roots. 

Besides the Jacobi-style version of EAI we may formulate the Gauss-Seidcl 
version of EAI, that is 



(i+i) (i) 

y) =y 



M(p{y?)) 



l-N{ P (yf))A 3 {y^ lV ^) 

■*3\y >y ) 2^ (i) (i+i) + (i) (i)- 
i=i Vj Ve e=j+i Vj 



The method, in the Jacobi version, is known to converge cubically for simple 
roots and linearly for multiple roots [40(. In the Gauss-Seidel version, con- 
vergence is slightly faster. In practice, good global convergence properties are 
observed; a theoretical analysis of global convergence, though, is still missing 
and constitutes an open problem. 

With the term vector iteration of the Ehrlich-Aberth method we refer to the 
step which provides the vector given the vector We use the term 

scalar iteration for indicating the single step performed on the generic scalar 
component of the vector yW. 

In the case of a scalar polynomial of degree N the cost of a scalar iteration 
is O(N) arithmetic operations. In this way, the cost of a vector iteration is at 
most 0(N 2 ) ops and is substantially reduced when most components have been 
numerically approximated, so that few scalar iterations must be performed in 
order to carry out the vector iteration. 



4 



The number of scalar iterations needed by the floating point implementation 
in order to find approximations which are exact roots of a slightly perturbed 
polynomial is, in practice, O(N) if the starting approximations are computed 
by means of the Newton polygon technique This technique is particularly 
effective when the polynomial has roots with moduli which are very unbalanced. 

Crucial aspects for an effective implementation of the EAI to matrix poly- 
nomials are 

1. the computation of the Newton correction p(x)/p'(x) given the value of x 
and of the input coefficients Pj, j = 0, . . . , k; 

2. a criterion for stopping the iterations; 

3. the choice of the initial approximations. 



2.1. Computing the Newton correction 

In the literature, methods based on some factorizations of P (x) were de- 
veloped to compute the Newton correction for functions that have the same 
zeros of p(x): e. g., t he method in [2l|, later proved to lack theoretical rigour 
and corrected in [19( • Other kinds of Newton-like approaches were presented in 



If one wishes to work with p{x) itself, a naive way to compute the Newton 
correction p(x)/p'(x) would be to evaluate first the coefficients of the polyno- 
mial p(x), say, by means of the evaluation- interpolation technique, and then to 
apply right after the Ehrlich-Aberth method to the scalar equation p(x) = 0. 
This approach would however come across numerical problems due to numerical 
instability and to overflow and underflow situations encountered in the compu- 
tation of determinants. 

It is therefore wise to conceive a strategy to avoid the explicit calculation of 
the coefficients of p(x). 

An effective way rests upon the well-known Jacobi's formula for the differ- 
ential of the determinant of any invertible square matrix A: 

d(ln(det A)) = tr(A~ 1 dA), (4) 

where tr denotes the trace. This way, we obtain the following expression for the 
derivative p'(x) 

= ddet(P(x)) = , dPp_ 

ax ax 
This formula allows us to evaluate the Newton correction p{x)/p'{x), which 
is the centerpiece for the EAI, without explicitly calculating p(x): 

P{X),AX) = tr (Pixy iP<(x)Y (5) 

An evaluation of P(x) and P'(x) by means of Horner's method, followed by 
a numerical matrix inversion, allows to compute the trace of in 
0(kn 2 + n?) operations. 
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We mention that, even though we independently formulated it, we found 
out later on that the possible use of §4§ in a numerical method for PEPs had 



been already suggested in [12|, |22|. However, in these instances it was proposed 
to use it to apply the Newton method to approximate each single eigenvalue in 
sequence, mentioning the possibility to use an implicit deflation of previously 
found roots [3(| , in order to avoid that the method converges twice to the same 
eigenvalue. This leads to a formula akin to ([2]), (|3]), with the difference that 
the summation in the term A is performed only up to the number of roots 
that have already been approximated. This is a crucial detail, because such a 
sequential implementation of the Newton method does not seem to achieve the 
same efficiency with respect to the EAI. 

2.2. Stopping criterion 

At the generic zth vector iteration it is crucial to decide whether the update 
of the jth component of the vector must be performed or the scalar 

iteration in that component must be halted. 

Observe that, if £ is a root of p(x), that is det P(£) = 0, then as the approx- 
imation x gets close to £ the matrix P(x) becomes ill-conditioned. This makes 
quite natural to stop the iterations if the reciprocal of the condition number 
fjt(P(x)) is less than a prescribed tolerance t\. This criterion makes sense if 
the eigenvalue that we want to approximate is semi-simple. In the case of a 
defective eigenvalue A with Jordan chains of length at most m, in view of the 
results in [39(, it is more convenient to stop the iterations if the reciprocal of 
fjt(P(x)) is less than r™. This latter condition is hard to implement since it is 
not easy to evaluate numerically the length of the Jordan chains of a matrix 
polynomial. Using the former stopping criterion may lead to a premature halt 
of the algorithm in the case of defective eigenvalues. 

As an alternative to the previous stopping condition, following |46j . define 

a ( x ) = X^=o \ x ^\- If Vj ^ s n °t an eigenvalue of P{x) then the quantity 

ri{yf) = (\\(p{yf)y 1 h{i + oc{yf))j 

measures the backward error for the approximation y?\ and can be cheaply 

evaluated during the EAI. The iteration can be halted when 77(2/^ ) is smaller 
than a given tolerance. 

For simple eigenvalues, no significant differences emerged between the two 
alternative possibilities. Therefore, our default choice was in favour of the cri- 
terion based on the condition number. 

It is also convenient to add, with the "or" logic operator, the following 
condition 

W{yf)l^-^(yf)My {i) ,y {i+1) )\ < -*\v?\ (6) 

where ti is a given tolerance. This condition says that the computed correction 
is too tiny and would not change the significant digits of the current approxi- 
mation. 
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2.3. Choosing initial approximations 

As pointed out in [l[ 0, [l5[ , practically effective choices of initial approxi- 
mations for the EAI are complex numbers equally displaced along circles. For 
instance, in [l[ it is proposed to choose initial approximations displaced along 
a circle centered at the origin of sufficiently large radius so that it contains all 
the roots. In [l5| the radius of the circle is suitably chosen. This strategy does 
not work effectively for polynomials having zeros with very large and with very 
small moduli. In |3| this drawback is overcome by considering different circles 
centered at the origin of suitable radii. The computation of these radii relies on 
the Rouche theorem. 

Here we try to extend this technique to a certain extent. We recall that, 
according to the Rouche theorem, if s(x) and q{x) are two polynomials such 
that 

\s(x)\ > \q{x)\, for |x| = r, 

then s(x) and s(x) + q(x) have the same number of roots in the open disk {z S 
C : |x| < r}. Applying this property with s(x) = x m and q(x) = p(x) — s{x), 
for < m < N : implies that if r m > X^jLo jjtm l a il r "' then the polynomial 
p{x) has m roots in the open disk of center and radius r. This property 
is at the basis of the criterion described in |3(, based on the Newton polygon 
construction, for choosing initial approximations equidistributed along different 
circles centered in 0. 

In order to extend this criterion to the case of matrix polynomials we need 
a generalisation of the Rouche theorem to matrix polynomials. We report the 
following result of [34| which we rephrase in a simpler way better suited for our 
problem. 

Theorem 1. Let S(x) and Q(x) be matrix polynomials and let r be a positive 
real. If S(x)* S(x) — Q(x)*Q(x) is positive definite for \x\ = r, then the polyno- 
mials det S{x) and det(S'(x) + Q(x)) have the same number of roots of modulus 
less than r. 

The following result is an immediate consequence of the above theorem ap- 
plied to the polynomial P(x) of ([T]) with S(x) = x m P m and Q(x) — J2i=o i^m x% Pi- 

Corollary 1. Assume that 

k k 

P^P m r 2m -( J2 E ^>0, for\x\=r, (7) 

j=0,j^m j=0,j^m 

where Ay- B means that A — B is positive definite. Then the matrix polynomial 
P(x) has mk eigenvalues in the open disk of center and radius r. 

Observe that if detP m = then condition (0 cannot be verified. In fact, 
the vector v such that P m v — would be such that 

k k 
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which is absurd. 

In particular, if det P k ^ the above corollary, applied with m = k, implies 
that all the eigenvalues of P(z) are included in the disk of center and radius 
r provided that 

k-1 k-1 

r 2fe P fc *P fc -(^P;^)(^P^>0, for|x| = r. (8) 

3=0 3=0 

Observe that the latter condition is implied by 

fe-i 

r 2k P* k P k y ^VP/P, +/^p(|P,-nP| + \Pim\). (9) 

3=0 3>i 

Similarly, applying Corollary [1] with m — provides a disk where P(x) has 
no eigenvalues. 

As an example of application, consider the 5x5 quadratic matrix polynomial 
P(x) = Ax 2 + Bx + A T , where B is the tridiagonal matrix defined by the entries 
[1, 2, 1], and A is the matrix with diagonal entries 100, 1, 1/1000, 1/100000, su- 
perdiagonal entries equal to 1 and with zero entries elsewhere. The eigenvalues 
of P(x) have moduli 2. 0050e+05, 1.4969e+03, 1.0000e+00, 1.0000e+00, 
1.0000e+00, 1.0000e+00, 6.6805e-04, 4 . 9874e-06. The criterion based on 
the above corollary in the form Q yields the bound 4.4e — 6 < |x| < 2.24e5 
which is quite good. Applying condition ([S]) yields the bounds 1.96e — 6 < \x\ < 
5.1e5 which is still good. 

Similar results have been obtained in [44| in the framework of tropical alge- 
bras. A different heuristic approach, which relies on Corollary [TJ is to select the 
values of the radii by considering the inequality ||P m ||2^ m > J2i=o ll-f'ilh''* 8 

in place of ([7]). This strategy, applied in the form r k > J2i=o \\Pi\\2 ri leads 
to the criterion based on computing the Newton polygon of the polynomial 
Si=o II -Pi II 2^*. This is the default choice of the starting approximations per- 
formed in our implementations. 

2.4- A posteriori error bounds 

In the case of a scalar polynomial p(x) of degree N, given a set of approxi- 
mations xi,. . . ,xn to the roots of p(x) it is possible to prove that [45| the set of 
disks Di — D{xi,ri) of center Xi and radius = n\p(xi)/(pN YljLi j*a( x i~ x j))\ 
is such that 

1. the union of the disks contains all the roots of p(x) 

2. each connected component formed by the union of, say, c overlapping 
disks, contains c roots of p(x). 

The set formed by Di, i — 1, . . . , N with the above properties is said set of 
inclusion disks. 

In the case of a matrix polynomial P(x) where det Pk ^ 0, it is quite cheap 
to compute a set of inclusion disks. In fact, if P(x) = ULU is the PLU fac- 
torization of P(x), then |p(x)| = |detP(x)| = Ilj=il u jjli wnere U = (uij). 
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Moreover, the leading coefficient pn of detP(x) coincides with detPfc which 
can be computed once for all. Observe that the LAPACK routine zgesv which 
solves a linear system with the matrix P(x), used to compute the Newton cor- 
rection l/trace(P(x) _1 P'(a:)), applied with x = Xi, provides at a negligible cost 
also the radius r». 

The availability of a set of inclusion disks enables one to perform a cluster 
analysis. In fact, once an isolated disk has been detected, we have isolated 
a single eigenvalue of the matrix polynomial P(x). Once we have detected a 
set of c overlapping disks isolated from the remaining inclusion disks, we have 
detected a cluster formed by c eigenvalues of P(x). 

A different a posteriori error bound can be obtained by using a classical 



result [17| . The disk of center Xi and radius f< = n\p(xi)/p'(xi)\ contains a root 



of the polynomial p{x). However, the set of disks obtained in this way does not 
fulfill properties 1 and 2 of the set of inclusion disks. It is worth pointing out that 
the computation of f, is inexpensive since the Newton correction p(xi)/p'(xi) is 
computed by the EAI. Moreover, this a posteriori error bound still holds if the 
leading coefficient Pfc is singular. 



2.5. Multiple eigenvalues 

Computational difficulties may be encountered in the case of multiple eigen- 
values. In fact, the rate of convergence for multiple eigenvalues is linear, with 
respect to the cubic behaviour for simple eigenvalues. Moreover, for defective 
eigenvalues the standard stop condition may lead to a premature halt. For this 
reason, if it is possible to detect a priori multiple eigenvalues, it is advisable to 
deflate them; if it is not possible to spot all of them theoretically, even lower 
bounds on the multiplicity are very helpful. If multiple eigenvalues are not pre- 
dicted theoretically, one must rely on the cluster analysis to identify them and 
modify accordingly the stopping criterion. 

A common situation that leads to multiple eigenvalues is met when the 
extremal coefficients are rank-deficient matrices. In this case, and/or oo have 
multiplicity greater than or equal to 1. This situation can be circumvented to 
a certain extent. 

In the case of to eigenvalues at infinity, one may just start with an approxi- 
mation vector y of length nk — to, acknowledging that the determinant p{x) has 
in fact degree nk — to; if there are to zero eigenvalues it is possible to set to zero 
to components of the vector y(°> avoiding to update them. 

The number of null singular values of Pq provides a lower bound to the 
number of null eigenvalues of P{x). Similarly, the number of zero singular 
values of Pfc provides a bound to the number of eigenvalues at infinity. This way, 
the precomputation of the SVD of Pq and P^ may increase the performance of 
the EAI. Equivalently, one may perform any rank- revealing factorization (e.g., 
QR) instead of the SVD. Sometimes, the structure of the coefficients allows 
to achieve better bounds (e.g., if the same rows/columns of many consecutive 
extremal coefficients are zero). 

In our implementation, the rank of the extremal coefficients is tested. If it is 
less than n, it is also checked if Pq and Pi (resp., Pk and Pc-i) share any common 
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zero row/column. Thus, any manifest presence of zero and infinite eigenvalues 
is exploited, forcing deflation of all the guaranteed roots. Moreover, if the test 
detects the presence of eigenvalues at (resp., oo), in order to avoid a premature 
stop for other undetected eigenvalues at (resp., oo) the stopping criterion is 
made stricter. The stronger stop condition requires that, for eigenvalues smaller 
(resp., larger) than a given bound, either the relative correction criterion ([B]) is 
satisfied with tolerance or the relative correction criterion (j6)) is satisfied with 

1/2 

tolerance r 2 and, simultaneously, the reciprocal condition number criterion 
is satisfied with tolerance t\. This heuristic device worked very effectively in 
our experiments, leading to satisfying results also in problems with multiple 
eigenvalues at either zero or infinity (see Section 01 . 

If the leading coefficient Pk is singular and if the degree of p(x) = det P(x) 
is not available together with the leading coefficients of p(x), then it is not 
possible to generate a set of inclusion disks and to perform a cluster analysis. 
However, in this case we may apply an effective technique based on a rational 
transformation of the variable x. For instance, the variable x is replaced by the 
Mobius function x = x(z) = (az + f5)/(^z + S) such that aS — 7/3 7^ 0, and the 
polynomial P(x) is replaced by the polynomial Q(z) = (72 + S) k P(x(z)). This 
way the eigenvalues at infinity of P(x) are mapped into eigenvalues of Q(z) at 
—S/j. Moreover, Q(z) has no eigenvalues at infinity provided that a/7 is not 
eigenvalue of P{x). The substitution of variable can be performed implicitly 
without actually computin g th e coefficients of Q(z) except for Qk- We refer the 
reader to Section [3] and to [37J for more details. 




2. 6. Linearization as a possibility 

In the paper [l3| the computation of the Newton correction was carried out 
by first linearizing the polynomial by means of companion-like pencil, and then 
by evaluating the trace by means of an LQ factorization. This approach has 
a computational complexity of 0(n 3 + kn 2 ) ops per scalar iteration, which is 
the same that is achieved by our algorithm. However, this cost can be reduced 
to 0(n 2 k 2 ) ops if the matrix pencil is reduced to triangular- Hessenberg form 
before the Jacobi formula is applied. This fact is the main advantage of using a 
linearization. However, linearization techniques, if not properly used, may lead 
to an undesired increasing of the eigenvalue condition numbers 1S\. 

3. The case of structured polynomials 

The EAI is particularly suited to deal with matrix polynomials endowed 
with specific structures of the matrix coefficients. We are interested in matrix 
structures which induce particular symmetries on the location of the eigenvalues. 
Polynomials of this kind are encountered in the applications and include, for 
instance, palindromic, T-palindromic, symplectic and Hamiltonian polynomials. 

Customary PEP-solving algorithms, such as the application of the QZ to 
any suitable linearization of the polynomial, are not able to fully catch these 
symmetries of the spectrum. In the literature, there are specific matrix methods 
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that achieve this goal. The EAI enables to exploit the additional information 
both in the computation of the Newton correction and in the choice and in the 
management of the (initial) approximation of the roots in view of the structure- 
induced symmetries. We will see this later on. 

Assume that the structured PEP is such that the eigenvalues appear in pairs 
{x, f{x)}, with f(f(x)) — x Vx. A naive adaptation of the EAI to this property 
would be to apply ([2]) or ((3j) updating only the first half of the components of the 
vector y and simultaneously imposing j/W = /(i/ l ~ rlfe / 2 )), i = nk/2 + 1, . . . , nk. 
In numerical experiments, this approach seems not to be always efficient in 
terms of number of scalar iterations needed for numerical convergence. This 
motivates the design of more sophisticated structured variants of the EAI, that 
we are going to describe in the following. 

Before analyzing the various classes of structured matrix polynomials, we 
recall some basic definitions of special matrices. 

An n x n square matrix A is said to be symmetric if A T = A and skew- 
symmetric if A T — ~A. Let n — 2m. The matrix A is said to be Hamiltonian 
if it is such that A T J — — J A where J is the matrix ( _°j ^g*); A is said to be 
skew-Hamiltonian if it is such that A T J = J A; A is said to be symplectic if it 
is such that A T JA = J. 

Every skew-Hamiltonian matrix can be obtained as the square of a Hamil- 
tonian matrix, and conversely the square of a Hamiltonian matrix is always 
skew-Hamiltonian [llj. Symplectic matrices are exponential of Hamiltonian ma- 
trices. 



3.1. Skew-Hamiltonian and even- dimensional skew- symmetric 

A skew-symmetric polynomial P{x) is a polynomial whose coefficients Pj, 
for j = 0, . . . , k, are skew-symmetric constant matrices. If the coefficients have 
even size n = 2m, we say that P{x) is an even-dimensional skew-symmetric 
polynomial. A skew-Hamiltonian polynomial is defined as a polynomial whose 
coefficients Pj are all skew-Hamiltonian matrices. Classical eigenvalue problems 
for skew-Hamiltonian matrices [48[ are a special case of skew-Hamiltonian PEPs. 

These two classes of polynomials are closely related, because multiplication 
by J maps one class onto the other. A common feature is that the spectrum of 
any polynomial in these two classes contains only eigenvalues of even multiplic- 
ity. In fact, the determinant of a matrix polynomial P(x) belonging to these 
two classes can be written as 



p(x) = detP(x) = q(x) ■ q{x), 

for a suitable polynomial q(x). For the special case of a real skew-symmetric 
matrix pencil, a proof was given in [24[ where a special Kronccker form was 
derived. The more general case comes from classical results on determinants 



35| . Let us give here a simple proof of the statement for an even-dimensional 



skew-symmetric complex matrix polynomial using modern terminology. 
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Proposition 1. Let P(x) = —P(x) T be a 2m x 2m skew- symmetric matrix 
polynomial. Then p{x) — detP(x) = q{x) ■ q(x) for some scalar polynomial 
q(x). 

Proof. We shall prove the proposition by induction on m. For m = 1 the 
statement is obvious. Suppose now that any (2m — 2) x (2m — 2) skew-symmetric 
polynomial has the desired property. Let II be a 2m x 2m permutation matrix 
and let Q(x) :— ILP(x)II T . Suppose that II is such that 



Qo(x) :=Q(1:2,1:2) =: 



r(x) 
-r(x) 



is nonsingular, where r{x) is a suitable nonzero scalar polynomial. Notice that 
such an assumption can be safely made because if that was false for any II 
then P(x) = so p(x) = and there would be nothing to prove. Now let 

Q(x) = (y^2(x) T Qi^Or))' wnere tne polynomial matrices A(x) and Q\{x) have, 
respectively, dimensions 2 x (2m — 2) and (2m — 2) x (2m — 2); also, let p{x) := 
r{x) m - 1 . Define the rational function S{x) := Qx{x) + ^A(x) T ( \ J) A(x). 
Clearly, r(x)S(x) is a (2m — 2) x (2m — 2) skew-symmetric matrix polyno- 
mial; therefore, by the inductive hypothesis, det S(x) — where 9{x) is a 



suitable scalar polynomial. Moreover, S(x) is the Schur complement of Qo(x) 

p(x) 

r(x) 



Thus, p(x) = ^~7^p^ i so p(x) is the square of some scalar rational function 
q(x) = f/^)L ■ Since p(x) is a polynomial, q(x) must be a polynomial as well. □ 



This is a particularly useful property which can be fully exploited by the 
Ehrlich-Aberth method. In fact, instead of applying the EAI to the polynomial 
p(x) of degree 2mfc, one can apply the EAI to the polynomial q(x) of degree mk 
even though q{x) is not explicitly known. 

More precisely, since p'(x)/p(x) — 2q'(x)/q{x), one can compute the Newton 
correction q(x)/q'(x) by means of 

q(x)/q'(x) = 2p(x)/p\x) = 2/tr(P(x)- 1 P' (x)). 

This way, the length of the vector of the approximations y in or in 
([3]) is reduced from 2mk to mk, moreover, the skew-Hamiltonian or the skew- 
symmetric structure of the coefficients can be exploited in the computation of 
P(x)- 1 P'(a;). 

3.2. Palindromic and symplectic 

The polynomial P{x) is called purely palindromic if RevP(a;) = P{x), where 
the reversal polynomial Rev(P(x)) is defined by RevP(x) :— x k P(x~ 1 ). The 
polynomial P(x) is called T-palindromic if RcvP(.t) = P(x) T . Both these 
structures induce a symmetry (x, 1 /x) in the spectrum. There is a vast literature 
on this kind of structure; see, e.g., ijj 25|, [27], and the references therein. 



The same structure appears in the standard eigenvalue problem for a symplectic 
matrix 101. 
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For this class of PEPs the change of variable z := x + 1/x is useful. In [13J it 
was shown that the use of a non-standard polynomial basis, called the Dickson 
basis, leads to a suitable linearization of the purely palindromic case. Moreover, 
it was shown that if P(x) is T-palindromic then it is possible to build a new 
skew-Hamiltonian matrix polynomial M{z) such that det M{z(x)) = p(x)-p(x); 
the Dickson basis was then used to obtain a useful linearization. In the follow- 
ing, we will show how to avoid the explicit use of the Dickson transformation. 
This has the advantage that there is no potential loss of accuracy for very large 
eigenvalues unlike the case of the algorithm in [l3[ where the linearization in- 
troduces unwanted defective eigenvalues at infinity which may create numerical 
problems (although such problems can be effectively amended by a structured 
refinement [bfjjV 

If nk is even, then by means of simple formal manipulations one may show 
that q(z) := x{z)~ nk / 2 ■ p(x(z)) is a polynomial in z, where x(z) — (z + 
\/ z 2 — 4)/2 or x(z) = {z — V z 2 — 4)/2, i.e., z(x) is one of the two branches 
of the inverse function of z(x) — x + 1/x. 

Moreover, taking derivatives in the latter equation leads to an explicit ex- 
pression for the Newton correction q{z)/q'{z) given in terms of p(x)/p'(x) 

W) = vm^iwr p ' {x),p[x) = t ^ lp '^- ^ 

This equation enables one to apply the EAI to the polynomial q(z). Once 
its roots Zi, . . . , z nk / 2 have been computed, the eigenvalues of P(x) are given by 
the pairs (xi, 1/xi) which are the roots of the quadratic polynomial x 2 — ZiX+ 1. 
This approach has the advantage to work with an approximation vector of half 
the size and to deliver the solution as pairs (x, 1/x). 

It is important to point out that the applications z — > x — (z ± V 1 z 2 — 4)/2 
is ill conditioned at z = ±2. Therefore, loss of accuracy is expected near x = 
±1. In this refinement step is advisable. Such a refinement may be 

implemented by an unstructured version of the EAI, by the naive strucutred 
EAI, or with other structured refinement methods (14] . 

If nk is odd, then —1 is necessarily an eigenvalue of the palindromic PEP 
and there is no need to approximate it. To calculate approximations of the 
remaining nk — 1 eigenvalues, there are two possible strategies. 

As a first possibility, one may consider the new matrix polynomial Q(x) = 
(x+l)P(x) which has even degree. The eigenvalues of Q(x) are those of P(x) and 
the eigenvalue —1 with multiplicity increased by n. Therefore, the previously 
described technique can be applied. Only nk — 1 roots of det(Q(x)) are needed, 
because n + 1 roots are a priori known to be equal to —1. Thus, one could apply 
the EAI © or ^ with an approximation vector y of n(k + 1) components of 
which n+1 are set equal to —1 in order to immediately achieve implicit deflation 
of the roots; or, working in the variable z in order to extract the structure, the 
SEAI (fTU|) can be used setting (n + l)/2 starting points equal to —2. 

A second possibility is to set x :— w 2 and to consider the eigenvalues of 
the polynomial Q(w) = P(x(w)). The scalar polynomial q(w) := det Q(w) has 
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2nk roots, which are the square roots of the solutions of the original equation 
p(x) = that we have to solve. In particular only 2nk — 2 roots are to be 
determined, since q(w) = has two known solution at w — ±i. It is useful to 
set z := (w + 1/w) 2 = x + l/x + 2. Defining 

~(w) ■= gM 

it is easy to check that r(z) := q(w(z)) is a polynomial in z. Therefore we 
may restrict the attention to computing the roots of r(z). Once they have 
been computed, the evaluation of the function w(z) at these roots provides the 
roots of q(w) . The evaluation of x(w) at these latter roots yields the sought 
eigenvalues of P(x). In order to compute the roots of r(z) we may apply the 
EAI to the polynomial r(z). The following equations provides a tool to compute 
the Newton correction r(z)/r'(z) needed by the EAI. 

r(z) _ 2w(l - 1/w 4 ) 

r'(z) q'(w)/q(w) — [(nk + l)w 2 + nk — l]/(u> 3 + w) ' 

or in terms of the original variable x 

r(z) _ 1 - l/x 2 

r'(z) ~ p'(x) /p(x) - [(nk + l)x + nk - l]/(2x 2 + 2x) ' 

At the moment we have no clear elements to say which of the two possibilities 
is more convenient. We plan to investigate in this direction. 

We conclude this subsection mentioning that also antipalindromic and anti- 
T -palindromic polynomials (RevP(a;) = —P(x) and RevP(x) = — P(x) T ) have 
a {x, l/x} symmetry. Their determinants are pure palindromic if n is even 
and antipalindromic if n is odd [26]. The former case is exactly the same as 
above. The latter case is also easy, because a scalar antipalindromic polynomial 
is always equal to x — 1 times a scalar pure palindromic polynomial. Moreover, 
it is possible to prove [26| that 1 is always a root of a scalar antipalindromic 
polynomial, and —1 is always a root of even-grade antipalindromic polynomial, 
so according to the grade there arc cither one or two exceptional eigenvalues 
with odd multiplicity. Therefore, it is easy to extend our technique to this class. 



3.3. Hamiltonian/skew-Hamiltonian and even/odd 

An even (odd) polynomial is such that Pj is symmetric for all even (odd) 
values of j and is skew-symmetric for all odd (even) j. Similarly, the coefficients 
of a Hamiltonian/skew-Hamiltonian polynomial are, alternatively, Hamiltonian 
and skcw-Hamiltonian matrices. The classes of even-dimensional even/odd poly- 
nomials are easily mapped onto the classes of Hamiltonian/skew Hamiltonian 
polynomials by a multiplication by J. Amongst the huge literature on these 
classes of polynomials see, for instance, 31 , H,[33|,|43| and the references therein. 



Classical eigenvalue problems for Hamiltonian matrices [48| are a special case 
of skew-Hamiltonian PEPs. 
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The matrix polynomials belonging to these classes have eigenvalues coming 
in pairs (x, —x). In particular, if nk is odd, then either x = (if Po is skew- 
symmetric) or x — oo (if Pk is skew-symmetric) is an eigenvalue. Notice that 
this is never the case for Hamiltonian/skew-Hamiltonian polynomials, because 
they are only defined for even n. 

Let z := x 2 . Just like the T-palindromic case, also for even/odd polynomials 
it is possible to follow the ideas exposed in [l3| and build a new matrix polyno- 
mial M(z) whose determinant is equal to p(x(z)) -p(x(z)). The following result 
demonstrates the way it can be done for an even polynomial. 



Proposition 2. Let P(x) be an even matrix polynomial, and let z — x 2 . Define 

B{z) f(.M)tf T (.(:)) and C{z) _^ P(x(.))-P J (x(.)) j fi0 that p {x) = 

B(x 2 ) + xC(x 2 ). Then M{z) := ^ *B(»? ) * s a skew-Hamiltonian matrix 
polynomial such that det M{z) = [p(x(z))} 2 . If ^ xq in C is an eigenvalue for 
P(x) associated with a canonical set of Jordan chains of length £x,...,£j~ then 
x^ is an eigenvalue for M(z) and its Jordan structure is the union of the Jordan 
structures of P{x) at xq and at —xq. 
Moreover: 

1. concerning eigenvectors associated with any finite nonzero eigenvalue xo, 
P(xo)vo = and P(—xq)wo — if and only if Vi — (xqVq ,Vq) t and 
wi = (—xqWq,Wq) t are two linearly independent eigenvectors such that 
M{xl)vi = = M{xl)w\; 

2. M(z) is a In x In matrix polynomial of degree degM = [(k + l)/2]; 

3. writing M(z) — X)j=o M Mjz^ , the relation Mj — (^J^ 2l +1 P p 2 1 ^ holds for 
all < j < deg M, where Pj = if j < or j > k; 

4. if k is odd, M(z) has at least n (respectively, n + 1) eigenvalues at infinity 
if n is even (respectively, odd). 



The proof can be obtained by adapting the arguments used in 13| for T- 
palindromic polynomials to the even case; we skip the details here. 

Similar results can of course be obtained for odd and Hamiltonian/skew- 
Hamiltonian matrix polynomials. Applying the EAI to M(z) allows us to extract 
the structure. An alternative approach, that avoids possible issues about loss 
of accuracy for very large eigenvalues (this time M(z) has extra infinite eigen- 
values only if k is odd), is once again the implicit use of the squaring transfor- 
mation. Namely, if nk is even (which is always satisfied for Hamiltonian/skew- 
Hamiltonian polynomials), then defining z := x 2 one finds that q(z) := p(x{z)) is 
a polynomial for x(z) — *J~z or x(z) = —sfz. Thus, p' (x) / (2xp(x)) = q'(z)/q(z), 
so that the Newton correction for the polynomial q{z) is readily available 

q(z)/q'(z) = 2xp(x)/p'(x) = 2x/tr(P{x)- 1 P' (x)) 

and the Ehrlich-Aberth algorithm can be implicitly applied to the polynomial 
q{z) in order to compute its roots Z\, . . . , z n k/2- This way, the roots of p{x) are 
readily available in pairs as {^fzl, —^fzl). 
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In other situations, one eigenvalue is necessarily either (if P(x) is odd 
and nk is odd) or oo (if P(x) is even and nk is odd); thus, there is no need 
to approximate it. It may also happen that there is one uncoupled eigenvalue 
at and another one at oo (e.g. if P(x) is odd, n is odd and k is even). 
In the case of an extra eigenvalue at 0, to approximate the other eigenvalues 
one can notice that q(z) :— p(y/z)/y/z is a polynomial and that q'(z)/q(z) = 
(l/2x 2 )(xp'(x)/p(x) — 1). This yields the Newton correction for q(z) as 

q(z)/q'(z) = 2x/(p'(x)/p(x) - 1/x) = 2x/(tr(P(x)- 1 P'( 2 ;)) - 1/x), z = x\ 

which enables one to apply the EAI to q(z) by using an approximation vector 
of length (nk — l)/2. As in the palindromic case, there is also the alternative 
option to consider the polynomial xP(x) which is even (odd) if P(x) is odd 
(even). The new polynomial xP(x) has n additional eigenvalues at that are 
known and can therefore be immediately deflated. 

3.4- Unified approach to any structure 

More in general, let C* := C U {oo} and let / : C* — > C* be any self-inverse 
function, that is f(f(xj) = x Vie C*. An example is the subclass of rational 
functions f(x) — ax +° , which are self-inverse whenever a 2 + be ^ 0. If we 
additionally require / to be analytic, having such a form is not only a sufficient 
condition, but it is also necessary (unless f(x) = x) for / to be self-inverse. This 
follows from the fact that Mobius functions (i.e., rational functions of degree 1) 
are the only automorphisms of C* . 

Suppose that, because of some structure in the coefficients of P(x), all eigen- 
values come in pairs {A, /(A)}. Eigenvalues such that A = /(A) are called ex- 
ceptional, and are allowed to appear with any multiplicity. 

Such a possibility justifies the requirements that 

1. f(x) is analytic, so that either it is the identity function or it has a finite 
number of fixed points, and 

2. there is a way to identify which exceptional eigenvalues, if any, appear 
with odd multiplicity. 

In fact, exceptional eigenvalues can otherwise become a problem. For instance, 
consider a real matrix polynomial associated to the non-analytic function f(x) = 
x*: the method meets problems in this case. The reason is that, since all the 
real line is exceptional, there is no way to state a priori which exceptional 
(i.e. real) eigenvalues, if any, appear without being part of a complex conjugate 
eigencouple. 

If / is analytic, the implicit change of variable method that we have described 
for the special cases f(x) = —x and f(x) = 1/x can be generalised in the 
following way. Suppose first that a ^ 0, and define z(x) := g£-±»£ = xf(x). Let 
x(z) denote any of the two branches of the inverse function of z(x). Then if there 
are no eigenvalues with odd multiplicity one can see that q(z) := f^z^zfjW^ * s 
a polynomial; therefore, the EAI can be applied to q(z) and the eigenvalues can 
be found inverting the rational function z(x). Otherwise (e.g. if nk is odd) there 
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must be some exceptional eigenvalues that can be treated with techniques akin 
to those described for the special cases previously considered. If on the contrary 
a = 0, let z(x) := cx ^ b = x+f(x). Once again if all the eigenvalues come out in 

couples then q(z) := P ^ X S^} 2 ^ s a polynomial; otherwise one can simply deal with 

x \ z ) 

the exceptional eigenvalues with known odd multiplicity by using techniques 
analogous to those described in the previous subsections. Notice that the fixed 
points of f(x) may lead to computational problems, since they are double roots 
in the equation z(x) = Q. Refinements of some kind are advisable there. See 
also [14 1 . 

The explicit method may also be extended to the general case. This is the 
subject of a future research project. 



4. Numerical experiments 

We have performed extensive numerical experiments in order to check the 
efficiency and the accuracy of our implementation of the EAI. Further tests have 
been performed to confirm the ability of our method to exploit structures in the 
coefficients and to respect structures in the spectrum when approximating it. 

4-1- Efficiency 

The complexity of the proposed algorithm is of order tn 3 + tkn 2 , where t is 
the total number of times that a trace computation is needed before (vectorial) 
convergence. There is empirical evidence that t heavily depends on the choice of 
the initial approximation. For the case of scalar polynomials, the use of suitable 
strategies [3j,|j| leads to a linear dependence of t with respect to the total number 
of roots. 

Experiments we made with our implementation, with starting points deter- 
mined by the Newton polygonal, suggest that this is also the case of the EAI 
applied to a matrix polynomial. This means that the computational complexity 
of the EAI is 0(kn 4 + k 2 n 3 ), leading to great computational advantages for 
k 3> \fn. As noticed in 12.61 if on the contrary k < yfn other more focused 
implementation of the EAI are possible, with cubic efficiency in kn. 

In order to confirm such predictions, we have compared our implementa- 
tion of the EAI and Matlab^'s QZ implementation polyeig on random matrix 
polynomials of high degree and small dimension. The experiments have been 
performed on a machine with CPU Intel Xeon 2.80GHz and system Linux De- 
bian 6.02. For very large values of k, we did not actually run polyeig due to the 
very large forecast computation times, but we extrapolated the times from the 
other experiments; in fact, when doubling the value of k we can expect that the 
running time of the QZ algorithm grows approximately with a factor 8. Such 
extrapolated values are marked with a * in the following tables. 



1 Matlab is a registered trademark of The MathWorks, Inc. 
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Computation times for n = 5 


Computation times for n = 2 


k 


EAI 


polyeig 


k 


EAI 


polyeig 


20 


0.062 s 


0.010 s 


50 


0.018 s 


0.015 s 


40 


0.121 s 


0.057 s 


100 


0.044 s 


0.064 s 


80 


0.312 s 


0.370 s 


200 


0.111 s 


0.369 s 


160 


0.920 s 


4.39 s 


400 


0.360 s 


4.35 s 


320 


2.92 s 


44.0 s 


800 


1.29 s 


51.9 s 


640 


10.3 s 


398 s 


1600 


4.76 s 


437 s 


1280 


38.1 s 


O(50 min)* 


3200 


18.4 s 


O(50 min)* 


2560 


148 s 


0(7 hours)* 




5120 


575 s 


0(2 days)* 



The values in the tables above are in agreement with our prediction that the 
computation time should asymptotically grow as k 2 . Moreover, the experimen- 
tation also confirms that, for a given value of n, the ratio of the time needed by 
EAI with respect to the time needed by the QZ algorithm exhibits an asymptotic 
growth that is approximately linear in k. This effect is taken to the extreme in 
the case n — 5, k = 5120. Had we used Matlab's polyeig, it would have taken 
several days of computation time on our machine to solve such a problem. Our 
implementation of the EAI gave the approximated eigenvalues in less than 10 
minutes. 

4-2. Accuracy 

In order to test the accuracy of our implementation we used the Matlab 
toolbox NLEVP[2j. This toolbox has been recently proposed by its authors as 
an interesting set of benchmark problems that may be used as a standard test 
for new methods for nonlinear eigenvalue problems. It contains data coming 
from practical applications as well as model problems known to have peculiar 
properties. 

Amongst the many nonlinear eigenproblems contained in NLEVP, we have 
selected all the square polynomial eigenproblems with n < 25k 2 . We discarded 
PEPs with a larger ratio n/k 2 because they could be better dealt with by a 
different implentation of EAI, via a preliminary linearization. The test suite 
selected with this criterion consists of 29 problems plus the 2 problems butterfly 
and wiresawl that, being structured, will be treated in the next subsection. 

In all the parameter-dependent problems in the library the default values of 
the parameters were selected. All methods were directly applied to the original 
matrices as saved in the library, without preprocessing them with any scaling. 
Forward errors are evaluated by comparing the approximations with either theo- 
retically known values, when available, or values computed in variable precision 
arithmetic (VPA) with Matlab's symbolic toolbox. 

The graphs below are in logarithmic scale. Whenever the absolute error for 
a certain eigenvalue A appeared to be numerically zero, i.e. it was less than A 
times the machine epsilon e = 2~ 52 ~ 2.22e — 16, we formally set it equal to 
4p Only absolute errors for the finite eigenvalues are shown in the figures. 
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For the 3 problems with k > 3, the eigenvalue forward errors where computed 
for both the EAI and the QZ method (as implemented in polyeig). Absolute 
errors for our implementation of the EAI are marked with a red * symbol, while 
absolute errors for polyeig are marked with a blue + symbol. For this set of 
experiments, we picked starting points on the unit circle. In our experience the 
order of magnitude of the forward error is not significantly affected by the choice 
of the starting points, even though for some problems other choices led to slight 
improvements (not discussed here). 




Fig. 1. Forward absolute errors for the problem orr sommerfeld 




Fig. 2. Forward absolute errors for the problem plasma drift 
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Fig. 3. Forward absolute errors for the problem relative pose 5pt 



For the 26 problems with k — 2, three methods were compared by computing 
their forward errors with the same method as above: polyeig (blue + symbol), 
EAI (red * symbol) and the software quadeig by Hammarling, Munro and Tis- 
seur [16J , specifically designed for quadratic PEPs (black x symbol) . Although 
we did not alter the coefficients given as input to any method, for most problems 
quadeig have performed scaling by the default settings of its internal algorithm, 
that prescribe scaling under certain conditions; see }16| . 



Fig. 4. Forward absolute errors for the problems acoustic wave Id (left) and 

acoustic wave 2d (right) 



Fig. 5. Forward absolute errors for the problems bicycle (left) and bilby (right) 
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Fig. 6. Forward absolute errors for the problems cd player (left) and closed 

loop (right) 



Fig. 7. Forward absolute errors for the problems dirac (left) and gen hyper 2 

(right) 



Fig. 8. Forward absolute errors for the problems hospital (left) and 
intersection (right) 
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Fig. 9. Forward absolute errors for the problems metal strip (left) and mobile 

manipulator (right) 



Fig. 10. Forward absolute errors for the problems omnicaml (left) and 

omnicam2 (right) 



Fig. 11. Forward absolute errors for the problems power plant (left) and qepl 

(right) 
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Fig. 12. Forward absolute errors for the problems qep2 (left) and qepS (right) 




Fig. 13. Forward absolute errors for the problems relative pose 6pt (left) and 

signl (right) 




Fig. 14. Forward absolute errors for the problems sign2 (left) and sleeper 

(right) 



23 



Fig. 15. Forward absolute errors for the problems spring (left) and spring 

dashpot (right) 



Fig. 16. Forward absolute errors for the problems wing (left) and wiresaw2 

(right) 

As can be seen by the figures above, the approximations of the EAI are 
competitive, and often more accurate than the approximations of the QZ. In 
some cases, the improvement is remarkable. We report in the following table 
the maximal relative error and the average relative error for all the finite (i.e. 
neither numerically zero nor numerically infinite) eigenvalues of the 29 consid- 
ered problems, and for both the EAI and the QZ. The average relative error is 
defined as the geometric mean of all relative errors; numerically zero relative 
errors have been counted as relative errors equal to e/2. The values reported 
for the QZ for quadratic problems correspond to the algorithm, picked between 
polyeig and quadeig, that achieved the best performance in terms of average 
relative error for the given problem. As can be deduced by the above pictures 
and coherently with the results on backward errors presented in [16| , such best 
performance was achieved generally, but not always, by the latter algorithm. 
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Problem 


Rcl. errors, EAI 


Rel. errors, QZ 




Max. 


Avg. 


Max. 


Avg. 


acoustic wave Id 


1.0e-14 


2.1c-16 


l.le-14 


1.7c-15 


acoustic wave 2d 


6/2 


6/2 


3.2c-15 


7.4c-16 


bicycle 


1.0e-15 


4.0c-16 


7.6c-15 


l.le-15 


bilby 


2.4c-14 


3.5c-16 


5.1c-15 


1.8e-15 


cd player 


5.3c-16 


1.2c-16 


4.0c-14 


3.3c-16 


closed loop 


e/2 


6/2 


3.4c-16 


1.5c-16 


dirac 


4.1e-14 


5.9e-15 


1.9c-13 


2.9e-14 


gen hyper 2 


2.5c-14 


9.3e-16 


2.4c-15 


4.9e-16 


hospital 


2.7e-15 


1.6e-16 


2.0c-14 


1.4c-15 


intersection 


4.8 e-9 


4.5 e-13 


1.0 


2.4e-8 


metal strip 


6.3e-16 


1.7c-16 


2.3c-15 


6.7e-16 


mobile manipulator 


6/2 


6/2 


5.1c-16 


5.1e-16 


omnicaml 


9.1c-ll 


1.2c-12 


4.3c-9 


6.4c-13 


omnicam2 


3.9c-10 


2.3c-15 


4.0c-9 


1.2e-13 


orr sommcrfeld 


5.0e-12 


9.1c-16 


4.8c-5 


2.8c-9 


plasma drift 


3.4e-13 


5.1c-16 


1.3c-ll 


4.7e-14 


power plant 


8.3e-14 


l.le-15 


6.1c-ll 


1.9e-13 


qepl 


8.9e-16 


1.7c-16 


1.8c-15 


5.1e-16 


qep2 


5.8c-9 


5.3c-ll 


2.2c-16 


1.9e-16 


qep3 


3.9c-9 


1.0e-14 


8.0c-10 


3.2c-14 


relative pose 5pt 


2.9e-14 


6.8e-15 


1.6e-14 


6.3e-15 


relative pose 6pt 


7.5e-14 


1.9e-14 


1.2c-13 


1.4c-14 


signl 


3.8c-8 


l.le-10 


5.0c-8 


3.9e-10 


sign2 


4.5c-14 


2.8c-15 


3.1c-13 


1.3c- 14 


sleeper 


8.0e-16 


2.8c-16 


1.8c-15 


5.6e-16 


spring 


6/2 


6/2 


1.7c-15 


3.5e-16 


spring dashpot 


5.6e-15 


3.1c-16 


2.8c-13 


1.4e-14 


wing 


6/2 


6/2 


1.2c-15 


8.1e-16 


wiresaw2 


6/2 


6/2 


2.4c-15 


9.2c-16 



As the results above show, the EAI was generally able to improve the accu- 
racy of the approximations with respect to the QZ method. The only problems 
where the EAI achieved an average performance worse than the QZ are QEP2 
(loss accuracy on the multiple eigenvalue 1 with respect to quadeig; polyeig has 
problems as well) and gen hyper 2. 

The problems in NLEVP do not have high degree, so the condition k 2 >• n 
is not met. Therefore, in contrast with the high degree case, for those problems 
using the EAI as a primary algorithm docs not bring advantages in term of 
computation time; on the contrary the implementation discussed in this paper 
is slower than QZ if n » k 2 . Also, in these cases a linearization-based version 
of the EAI would be more efficient than the polynomial-based implementation 
discussed in the present paper. However, it is worth noticing that the numerical 
experiments showed that very often it improved the accuracy achieved by polyeig 
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and/or quadeig; in many cases, the EAI is able to compute correctly all the digits 
of the eigenvalue up to machine precision. This suggests that, when k 2 < n, 
it is possible to use the EAI as a refinement algorithm in order to improve the 
approximations obtained by the QZ. Using such values as starting points offers 
of course a very good choice, lowering the number of overall scalar iteration 
needed before convergence and therefore improving the efficiency of the EAI. 



4-3. Structured case 

Discussions and analyses of the behaviour of different structured versions of 



the EAI applied to structured PEPs have already appeared in [13[ for palin- 
dromic polynomials and in .'lii for even/odd polynomials. 

Here we further verify the reliability of our method with two tests. The first 
test relies on the NLEVP problems butterfly and wiresawl, which are even 0. 
The next figures compare four algorithms applied to these problems: polyeig^s 
QZ (blue -I- symbol), unstructured EAI (UEAI, red * sy mbol), the structured 
matrix method URV applied to an even linearization [431 ] (black x symbol) , and 
a structured version of the EAI relying on the change of variable z = x 2 method 
(SEAI, green o symbol). 



Fig. 17. Forward absolute errors for the problem butterfly 
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Fig. 18. Forward absolute errors for the problem wiresawl 



It is clear from the figures above that for the problem butterfly the EAI 
was more accurate than the two matrix methods, but structured methods did 
not improve much the accuracy of each unstructured counterpart. This sug- 
gests that for these matrix polynomials the unstructured condition numbers for 
the eigenvalues are not much different from the structured condition numbers; 
therefore, structured methods do not improve much the accuracy, even though 
they improve the efficiency. 

In the second test, a matrix polynomial W(x) with n = 2 and k = 10 was 
built in such a way that its eigenvalues appear in couples of the form {A, t^t}- 
In order to devise a problem not too easy to solve numerically, the determinant 
of the polynomial was designed to be Wilkinson-like: det(W(x)) = const. ■ 0(x) ■ 
0(jrzy), 6{x) — x ■ nj=2( x ~ i)- The next figure shows the absolute errors of 
the computed approximations with respect to the known exact eigenvalues for 
three methods: QZ (polyeig, blue + symbol), UEAI (red * symbol) and SEAI 
relying on the change of variable z — x (green o symbol). Numerically zero 
errors were formally set equal to e/2. 
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Fig. 19. Forward absolute errors for the structured problem det W(x) = 



The following table reports the relative errors of the three methods consid- 
ered above for all the eigenvalues but (all the three algorithms detected the 
zero eigenvalue with an absolute error smaller than the machine epsilon). It 
also reports the relative error for the Matlab's function eig (without scaling) 
applied to a suitable linearization, chosen according to the prescriptions of jl 81 ] . 
Notice that here all the nonzero eigenvalues have modulus > lj_so the suggested 
(near-to-optimal) linearization in the space ML according to [18| would in p rin- 
ciple be the pencil in DL corresponding to the ansatz vector e\ (see [28| for 
further details). Unluckily, since one eigenvalue is zero, that pencil is not a lin- 
earization at all [28j. Following the suggestions on conditioning of [l8j and the 
theory on linearizations of (28|, we have therefore taken the slightly perturbed 
vector ei + 2~ 23 eio as an ansatz vector. The factor 2 -23 has been heuristically 
chosen picking the integer a < 52 that minimises the average relative error 
when applying eig to the linearizations in DL associated to the ansatz vectors 
d + 2-«ei . 
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Eigenvalue 


R. e., polyeig 


R. e., eig 


R. e., UEAI 


R. e., SEAI 


-1 


4.9c-9 


e/2 


e/2 


e/2 


11/9 


4.5c-2 


4.8e-8 


3.5c-9 


1.4e-13 


5/4 


8.1c-2 


2.2e-7 


l.le-8 


5.7e-13 


9/7 


8.9e-2 


4.2e-7 


4.0e-9 


1.4e-12 


4/3 


l.le-3 


4.1e-7 


l.le-8 


4.6e-12 


7/5 


9.9c-2 


2.2e-7 


7.1c-ll 


5.1c-12 


3/2 


8.9c-2 


6.3 e-8 


3.8c-10 


2.5c-12 


5/3 


3.9e-2 


8.9e-9 


1.8c-10 


l.lc-12 


2 


1.0e-10 


4.9c-10 


8.8c-12 


5.1c-14 


2 


7.3c-4 


2.3e-12 


6.7e-15 


4.7e-14 


3 


4.3e-6 


6.3e-12 


2.2e-14 


6.3e-14 


3 


6.5e-8 


5.8e-12 


1.5e-13 


6.7e-14 


4 


9.6e-8 


1.5e-10 


1.7c-12 


2.0e-12 


5 


1.3e-7 


8.1e-10 


2.9c-12 


5.9c-12 


6 


3.2e-7 


2.2e-9 


6.1e-12 


1.5e-ll 


7 


3.1c-6 


3.3e-9 


1.5 c-11 


1.6e-ll 


8 


4.6e-6 


2.7e-9 


1.5e-ll 


5.5e-12 


9 


2.9e-6 


l.le-9 


1.7c-12 


2.5e-12 


10 


4.2e-6 


1.7e-10 


3.0e-12 


6.9e-13 


Average 


3.8e-5 


9.6e-10 


5.5e-12 


4.1c-13 



We may conclude that on this structured problem the EAI outperforms the 
QZ method for what concerns accuracy. Apparently polyeig struggles quite a bit 
here (which is coherent with the results of [18( ) , while the strategy of [18j works 
better, but still worse than the EAI. Although there are some approximations 
that do not benefit from the use of the structured version of the Ehrlich-Aberth 
algorithm, the SEAI has an overall advantage in accuracy over the UEAI, besides 
the obvious efficiency advantage. The computation time for the SEAI was about 
one third of the computation time for the UEAI. 

5. Conclusions and lines of future research 

We have proposed and tested a generalisation of the Ehrlich-Aberth method 
to polynomial eigenvalue problems. Both theoretical arguments and numerical 
experiments show that the Ehrlich-Aberth algorithm is more efficient than cus- 
tomary method for high-degree PEPs, since its complexity is only quadratic in 
k. Numerical experiments also suggest that in many situations the new method 
provides more accurate approximations, and therefore it may also be used as a 
refinement method for low-degree PEPs. 

On the other hand, problems arise in the treatment of multiple eigenvalues. 
Current research is focused on this issue. Currently, our algorithm exploits a 
heuristic device to deal with multiple eigenvalues at or oo, which is the most 
common occurrence in practice. For the NLEVP problems, very good results 
were obtained also when multiple eigenvalues were present. 
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Another line of current research regards the possible use of other root-finding 
algorithms. For instance, we can cite the modified Ehrlich-Aberth iteration 
1H , the Durand-Kerner iteration [29|, [Z(| or simultaneous root finders based 



on higher order methods in the Householder family, e.g. the Halley method [401 ] - 
It is advocated 2§, 4(| that some of the above mentioned method have order 



of convergence higher than 3. However, in practice we did not see a significant 
improvement on the total number of scalar iterations with respect to the EAI; 
sometimes, performances were definitely worse than the EAI. Further work is 
needed to compare the various possibilities in special cases like structured matrix 
polynomials. 

Important issues on which we plan to keep on working are improving the 
design of reliable stop conditions, a posteriori error bounds and choices of start- 
ing approximations. About the latter issue, in [|| some results of [44| on scalar 
polynomials will be generalised to matrix polynomials. 

Finally, we have proposed and tested a structured version of the algorithm 
that is able to catch certain structures in the spectrum. New numerical ex- 
periments on the matter confirm and extend the results of [l3|, [36| about the 
efficiency and the accuracy of the structured EAI. 



6. Acknowledgements 

The second author wishes to thank Federico Poloni who provided some help 
by delivering the eigenvalues computed by Matlab in VPA for some of the 
NLEVP problems. 



References 

[1] O. Aberth. Iteration methods for finding all zeros of a polynomial simulta- 
neously. Math. Comput, 27:339-344, 1973. 

[2] T. Betcke, N. J. Higham, V. Mehrmann, C. Schroder and F. Tis- 
seur. NLEVP: A Collection of Nonlinear Eigenvalue Problems. 
MIMS EPrint 2010.98, 2010. Available online on the webpage 
http: / /www.mims. manchester.ac.uk/research/numerical-analysis /nlevp.html| 

[3] D. A. Bini. Numerical computation of polynomial zeros by means of 
Aberths method. Numer. Algorithms, 13:179200, 1996. 

[4] D. A. Bini and G. Fiorentino. Design, analysis, and implementation of a 
multiprecision polynomial rootfinder. Numer. Algorithms, 23(2-3):127-173, 
2000. 

[5] D. A. Bini, L. Gemignani and F. Tisseur. The Ehrlich-Aberth method for 
the nonsymmetric tridiagonal eigenvalue problem. SIAM J. Matrix Anal. 
Appl, 27(1):153-175, 2005. 

[6] D. A. Bini, V. Noferini and M. Sharify. On some polynomial eigenvalue 
location theorems. In preparation. 



30 



[7] W. Borsch-Supan. A Posteriori Error Bounds for the Zeros of Polynomials. 
Numer. Math., 5:380-398 (1963). 

[8] R. P. Boyer and W. M. Y. Goh. Partition polynomials: asymptotics and 
zeros. Tapas in experimental mathematics, 99111, Contemp. Math., 457, 
Amer. Math. Soc., Providence, RI, 2008. 

[9] L. W. Ehrlich. A modified Newton method for polynomials, Comm. ACM, 
10(2):107-108, 1967. 

[10] H. Fafibcndcr. Symplectic methods for the symplectic eigenproblem. Kluwer, 
New York, 2002. 

[11] H. Fafibender, D. S. Mackey, N. Mackey and H. Xu. Hamiltonian Square 
Roots of Skew-Hamiltonian Matrices. Linear Algebra Appl., 287:125-159, 
1999. 

[12] W. Gander. Zeros of determinants of A-matrices. In Vadim Olshevsky and 
Eugene Tyrtyshnikov, editors, Matrix methods: theory, algorithms and ap- 
plications, Dedicated to the memory of G. Golub. World Scientific Publisher, 
2010. 

[13] L. Gemignani and V. Nofcrini. The Ehrlich-Aberth method for palindromic 
matrix polynomials represented in the Dickson basis. Linear Algebra Appl. 
2011, doi: 10. 1016/j.laa.2011. 10.035 

[14] L. Gemignani and V. Noferini. Modifications of Newton's method for even- 
grade palindromic polynomials and other twined polynomials. Submitted. 

[15] H. Guggenheimer. Initial approximations in Durand-Kerner's root finding 
method. BIT 26:537-539. 

[16] S. Hammarling, C. Munro and F. Tisseur. An Algorithm for the Com- 
plete Solution of Quadratic Eigenvalue Problems. MIMS EPrint 2011.86, 
available online. 

[17] P. Henrici. Applied and computational complex analysis. Vol. 1. Wiley- 
Interscience [John Wiley & Sons], New York-London-Sydney, 1974 

[18] N. J. Higham, D. S. Mackey and F. Tisseur. The conditioning of lineariza- 
tions of matrix polynomials. SIAM J. Matrix Anal. Appl, 28(4):1005-1028, 
2006. 

[19] N. K. Jain and K. Singhal. On Kublanovskayas approach to the solution 
of the generalized latent value problem for functional A-matrices. SIAM J. 
Matrix Anal. Appl, 20(5) T062-1070, 1983. 

[20] N. K. Jain, K. Singhal and K. Huscyin. On roots of functional lambda 
matrices. Comput. Meth. Appl. Mech. Engrg., 40:277-292, 1983. 



31 



[21] V. Kublanovskaya. On an approach to the solution of the generalized latent 
value problem for A-matrices. SIAM J. Matrix Anal. Appl., 7 :532-537, 
1970. 

[22] P. Lancaster. Lambda-matrices and vibrating structures. Pergamon Press, 
Oxford, 1966. 

[23] P. Lancaster, U. Prells and L. Rodman. Canonical structures for palin- 
dromic matrix polynomials. Oper. Matrices, 1:469-489, 2007. 

[24] P. Lancaster and L. Rodman. Canonical forms for symmetric/skew- 
symmetric real matrix pairs under strict equivalence and congruence. Lin- 
ear Algebra Appl, 406:1-76, 2005. 

[25] D. S. Mackey, N. Mackey, C. Mehl and V. Mehrmann. Numerical methods 
for palindromic eigenvalue problems: Computing the anti-triangular schur 
form. Numerical Linear Algebra with Appl, 16:63-86, 2009. 

[26] D. S. Mackey, N. Mackey, C. Mehl and V. Mehrmann. Smith forms of 
palindromic matrix polynomials. Electron. J. Linear Algebra, 22, 53-91, 
2011. 

[27] D. S. Mackey, N. Mackey, C. Mehl and V. Mehrmann. Structured polyno- 
mial eigenvalue problems: good vibrations from good linearizations. SIAM 
J. Matrix Anal. Appl, 28(4):1029-1051 (electronic), 2006. 

[28] D. S. Mackey, N. Mackey, C. Mehl and V. Mehrmann. Vector Spaces 
of Linearizations for Matrix Polynomials. SIAM J. Matrix Anal. Appl, 
28(4):971-1004 (electronic), 2006. 

[29] J. M. MacNamee. Numerical Methods for Roots of Polynomials. Studies in 
Computational Mathematics, 14, Elsevier (2007). 

[30] H. Mahley. Zur Iterativen Auflosung Algcbraisher Gleichungen. Z. Angew. 
Math. Physik 5:260-263, 1954. 

[31] C. Mehl. Condensed forms for skew-Hamiltonian/Hamiltonian pencils. 
SIAM J. Matrix Anal. Appl, 21:454-476, 1999. 

[32] V. Mehrmann and D. Watkins. Structure-preserving Methods for Com- 
puting Eigenpairs of Large Sparse Skew-Hamiltonian/Hamiltonian Pencils. 
SIAM J. Sci. Comp., 22:1905-1925, 2001. 

[33] V. Mehrmann and D. Watkins. Polynomial eigenvalue problems with 
Hamiltonian structure. Electr. Trans. Num. Anal, 13:106-113, 2002 

[34] Y. Mondcn and S. Arimoto. Generalized Rouche's theorem and its applica- 
tions to multivariate autoregressions. IEEE Trans. Acoust., Speech, Signal 
Processing, ASSP-28:733-738, 1980. 



32 



[35] T. Muir. A Treatise on the Theory of Determinants. Macmillan and Co., 
1882. 

[36] V. Nofcrini. The application of the Ehrlich-Aberth method to structured 
polynomial eigenvalue problems. Proc. Appl. Math. Mech., 11:919-922, 
2011. 

[37] V. Noferini. The behaviour of the complete eigenstructure of a polynomial 
matrix under a generic rational transformation. Submitted. 

[38] A.-W. M. Nourein. An iterative formula for the simultaneous determination 
of the zeros of a polynomial. J. Corn-put. Appl. Math., 1(4): 251-254, 1975. 

[39] N. Papathanasiou and P. Psarrakos. On condition numbers of polynomial 
eigenvalue problems. Appl. Math. Comput, 216(4):1194-1205, 2010. 

[40] M. Pctkovic. Point Estimation of Root Finding Methods. Lecture Notes in 
Mathematics, 1933. Springer- Verlag, Berlin Heidelberg, 2008. 

[41] B. Plestenjak. Numerical method for the tridiagonal hyperbolic quadratic 
eigenvalue problem. SIAM J. Matrix Anal. Appl., 28(4):1157-1172, 2006 

[42] A. Ruhe. Algorithms for the nonlinear eigenvalue problem. SIAM J. Numer. 
Anal, 10:674-689, 1973. 

[43] C. Schroder. Palindromic and Even Eigenvalue Problems - Analysis and Nu- 
merical Methods. PhD thesis, Technischen Universitat Berlin, April 2008. 

[44] M. Sharify. Scaling Algorithms and Tropical Methods in Numerical Ma- 
trix Analysis: Application to the Optimal Assignment Problem and to the 
Accurate Computation of Eigenvalues. Ph.D. Thesis, Ecole Polytechnique 
ParisTech, Specialite: Mathematiques Appliquees. 

[45] B. T. Smith. Error bounds for zeros of a polynomial based upon Ger- 
schgorin's theorems. J ACM, 17:661-674, 1970. 

[46] F. Tisseur. Backward Error and Condition of Polynomial Eigenvalue Prob- 
lems. Linear Algebra Appl., 309(l-3):339-361, 2000. 

[47] F. Tisseur, and K. Meerbergen. The quadratic eigenvalue problem. SIAM 
Rev., 43(2):235-286, 2001. 

[48] C. Van Loan. A Symplcctic Method for Approximating all the Eigenvalues 
of a Hamiltonian Matrix. Linear Algebra Appl., 61:233-251, 1984. 

[49] G. T. Wilson, The factorization of matricial spectral densities. SIAM J. 
Appl. Math., 23:420-426, 1972. 



33 



